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Abstract. We present benchmark problems and solutions for the continuum radiative transfer (RT) in a 2D disk configuration. 
The reliability of three Monte-Carlo and two grid-based codes is tested by comparing their results for a set of well-defined cases 
which differ for optical depth and viewing angle. For all the configurations, the overall shape of the resulting temperature and 
spectral energy distribution is well reproduced. The solutions we provide can be used for the verification of other RT codes. We 
also point out the advantages and disadvantages of the various numerical techniques applied to solve the RT problem. 

Key words, radiative transfer - stars: circumstellar matter - methods: numerical 



1. Introduction 

Observations show that many astrophysical sources such as 
young stellar objects (YSOs), post-AGB stars and active galac- 
tic nuclei (AGN) are surrounded by dust. Dust grains scatter, 
absorb and re-emit radiation originating from the primary en- 
ergy sources, thus modifying their spectral energy distributions 
(SEDs). Moreover many embedded objects cannot be directly 
studied in the visible, since dust may entirely obscure them at 
optical wavelengths. Their structure can be only inferred from 
the thermal dust emission. Therefore, modelling of their inten- 
sity and polarization maps as well as their SEDs is necessary. 
This can only be done by solving the radiative transfer (RT) 
equation (e.g. Yorke 1985). Analytical solutions for this equa- 
tion do exist only for the simplest cases, far from representing 
the complexity of dust-enshrouded objects. Hence, the devel- 
opment of sophisticated numerical RT codes is unavoidable. 

Early attempts for spherically symmetric configurations 
were performed by Hummer & Rybicki ( 119711 1. Scoville & 
Kwan ( 1976 1 and Leung ( 1976 1 including rough assumptions 
such as grey opacity and/or neglecting scattering. The first 
formal solution for the dust continuum in spherical geome- 
try was obtained by Rowan-Robinson (1980), directly inte- 
grating the RT equation, an operation known as ray- tracing. 
Since then many other codes treating 1-D slab or spherical 
configurations (e.g. Yorke & Shustov 119811 Lefevre et al. 



[19*521 Martin et al. "19831 Rogers & Martin [19*531 Henning 
1985; Groenewegen 1993, Winters et al. 1994) or the inverse 
RT problem (Steinacker et al. 2002b) have been developed. 
Nevertheless, it became soon clear that a 1-D geometry is often 
too restrictive. Distinct non-spherical features such as bipolar 
outflows (e.g. Bachiller 1996 ), bipolar reflection nebulae (e.g. 
Lenzen [T98""t and disks (e.g. McCaughrean & O'Dell 1996) 
are typical of many astronomical objects. Nowadays, multidi- 
mensional codes implementing different methods and numeri- 
cal schemes are being applied to treat the RT in such configu- 
rations. 

In contrast to hydrodynamical simulations, benchmark tests 
for radiative transfer computations are rare. The only prac- 
tical approach to test the reliability of RT calculations is to 
compare solutions of well-defined problems by several inde- 
pendent codes. This has been done for the ID case by Ivezic 
et al. dl997> . A benchmark project for 1-D plane-parallel RT 
and vertical structure calculations for irradiated passive disks 
is available on the web '. As for 1-D RT in molecular lines, a 
comparison of results from different codes has been performed 
by van Zadelhoff et al. (2002) 1 . Going from spherical symme- 
try to 2D and 3D spatial configurations, we add two or three 
more variables to the RT problem. Numerically, this implies 
10 4 or 10 6 more numbers to store when a decent resolution 



Send offprint requests to: I. Pascucci, pascucci@mpia-hd.mpg.de 



1 http://www.mpa-garching.mpg.de/PUBLICATIONS/DATA/radtrans/ 

benchmarks/ 

2 see also: http://www.strw.leidenuniv.nl/~radtrans/ 



2 



I. Pascucci et al.: The 2D Continuum Radiative Transfer Problem 



of 100 points in each variable is used. In addition, the geom- 
etry makes the solution of the integro-differential RT equation 
more complex. This explains why benchmark tests for 2D and 
3D configurations are lacking. It also implies that reaching an 
agreement to the level of ID RT computations using state-of- 
the-art computer equipment is unrealistic. A previous attempt 
to test 2D RT calculations has been made by Men'shchikov 
& Henning (1997). They compare results from their approxi- 
mate method with those of a fully-2D program (Efstathiou & 
Rowan-Robinson 1990 1 applying the same geometry. 

Here, we propose to test the behaviour of five different RT 
codes in a well defined 2D configuration, point out advantages 
and disadvantages of the various techniques applied to solve the 
RT problem and provide benchmark solutions for the verifica- 
tion of continuum RT codes. As modelling sources with high 
optical depth and strong scattering is the challenge of multi- 
dimensional RT codes, we explicitly include a test case at the 
limit of the current computational capabilities. In Sect. [2] we 
briefly introduce the RT problem and we define our test case. 
Sect. [3]is devoted to the description of the codes we used and 
to explain their differences. Solutions for the dust temperature 
and emerging SEDs are presented in Sect.@] In the last section 
we discuss our results. 

2. Benchmark Problem 

2.1. The RT problem 

Solving the RT problem means to determine the intensity 
I A (x, n) of the radiation field at each point x and direction n of 
the model geometry and at each wavelength A. This is achieved 
by solving the stationary transfer equation 

nV x h(x,n) = - [ K ab \A,x) +k™(A,x)] I A (x,n) 
+K abs (A,x) B A [T(x)] 

A- sca C/i x) r 

+ - d£l'p(A,n,n')I A (x,ri) 

\n J 

n 

+E A (x,n) (1) 

where K abs (/l, x) and K sca (/l, x) are the absorption and scat- 
tering coefficients of the particles, respectively. The quantity 
p(A, n, n') denotes the probability that radiation is scattered 
from the direction n' into n, Q is the solid angle, B A is the 
Planck function, and T is the temperature. The index A denotes 
that the quantity is defined per wavelength interval. E A (x, n) 
represents all internal radiation sources such as viscous heating 
or cosmic rays. For the sake of simplicity, we only consider one 
dust component of specific size and chemical composition. In 
addition, we do not discuss the polarization state of the radia- 
tion field and consider the intensity only. 

If spherical symmetry in the particle distribution and the 
sources of radiation is assumed, the integro-differential equa- 
tion Q becomes a function of 3 variables, already difficult to 
solve even for a given dust temperature T(x). In the case of spa- 
tial 2D configurations (axial-symmetric disks, tori), we have to 
deal with 5 variables. Moreover, the coupling between the radi- 
ation field and the dust temperature requires the simultaneous 



consideration of the balance equation for the local energy den- 
sity at point x 

oo oo 

J dA Qf s B A [T tad (x)] = J dAQ^j- J dQ! I A (x,n') (2) 

o on 

to calculate intensity and temperature self-consistently. Here, 
2 abs (/l) is the absorption efficiency factor, while T rad is the 
temperature arising from radiative heating. Additional heating 
sources can contribute to the temperature with 

T(x) = T md (x) + T heat (x). (3) 
2.2. Model definition 

We consider the general astrophysical case of a star embedded 
in a circumstellar disk with an inner cavity free of dust. We 
assume that the star is point-like, located at the center of the 
configuration and radiating as a black body at the same tem- 
perature as the Sun. The disk is made of spherical astronomi- 
cal silicate grains, having a radius of 0.12 pm and a density of 
3.6 gcrrT 3 (optical data are taken from Draine & Lee 1984', 
see also Fig. [Q. The disk radially extends to a maximum dis- 
tance of 1000 AU from the central star. Since the correct deter- 
mination of the sublimation radius is quite a difficult problem, 
we fix the inner radius to 1 AU. This guarantees a maximum 
dust temperature less than 1000 K, even in the case of high opti- 
cal depth. The density structure is that of a massless (in relation 
to the central star) Keplerian disk having no cutoff at a certain 
opening angle. This implies that the radiative transfer has to be 
simulated both in the optically thick disk and in the optically 
thin envelope. The disk geometry and density structure are sim- 
ilar to those described by Chiang & Goldreich J1997I 1999 1 and 
successfully applied to study passive disks around T Tauri stars 
(Natta et al. |2000l. The density distribution provides a steep- 
density gradient in the inner part of the disk which could give 
rise to numerical problems when solving the RT equation. This 
turns out to be an advantage for RT comparison since it allows 
to test the codes' behaviour under extreme conditions. The den- 
sity distribution we adopt has the following form 

p(r, z) = pa x /] (r) x f 2 (z/h(r)) (4) 

Mr) = (r/r d y u) 

f 2 (r) = exp(-7r/4x(z/h(r)) 2 ) 

Mr) = z d x (r/r d ) 1125 

with r being the distance from the central star in the disk mid- 
plane ( -\/x 2 + y 2 ) and z the distance from the midplane. Here r d 
is half of the disk outer radius (R ollt /2) and z d one fourth of r d 
(Rout/8). Note that the disk is slightly flared, i.e., the disk open- 
ing angle h{r)/r is exponentially increasing with the distance 
from the star. The term f\ provides the radial dependence of the 
density distribution. In protoplanetary disks, the volume den- 
sity is usually proportional to r~ a with a in the range (1.8H-2.8) 
(e.g. Wood et al.|2002 and Cotera et al. l2001> . For this bench- 
mark we use a = 1 in order to save CPU time. Both f\ and f 2 

3 downloadable from: 
http://www.mpia.de/PSF/PSFpages/RT/benchmark.html 
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Table 1. Model parameters 



Symbol 


meaning 


value 


M, 


Stellar mass 


1M 


o 

K, 


Stellar radius 


1 R 

1 Kg 


T, 


Stellar effective temperature 


5800 K 


^out 


Outer disk radius 


1000 AU 


Rm 


Inner disk radius 


1 AU 


Zd 


Disk height 


125 AU 


a 


Grain radius 


0.12yum 


Pg 


Grain density 


3.6 gcirT 3 


r v 


Optical depth at 550 nm 


0.1, 1, 10, 100 



remain unchanged while po is chosen so to define different opti- 
cal depths. We perform calculations for four values of visual (A 
= 550 nm) optical depth, namely t v =0.1, 1, 10, 100. The op- 
tical depth, as seen from the centre, is calculated along the disk 
midplane. Since most of the dust is confined in the midplane, 
the optical depths we refer to are the highest in each model. The 
test case t v = 100 is at the limit of our current computational 
capabilities. The resulting total dust mass for the model with 
r v = 1(100) is l.lxlO- 6 M (l.lxlO- 4 M ). The density struc- 
ture perpendicular to the disk midplane is shown for the same 
model in Fig. [5] The RT is calculated for 61 wavelengths be- 
ing distributed nearly equidistantly on a logarithmic scale from 
0. 12-2000 fim. These 61 wavelengths define the frequency res- 
olution of our computations. In Sect. 14.41 we also compare two 
Monte Carlo (MC) codes on a grid with two times more wave- 
lengths and we discuss the effect of the frequency resolution on 
the 2D benchmark. Since anisotropic scattering is not included 
in all codes, we consider the scattering as isotropic. Symbols 
and values of the model parameters are summarized in Tabled 
for more clarity. 

3. RT Simulations 

3.1. Methods 

Similar to hydrodynamical simulations, we can distinguish par- 
ticle (Monte Carlo) and grid-based methods to solve the RT 
equation numerically (Henning 2001 1. 

In MC simulations the radiation field is partitioned in 
equal-energy, monochromatic "photon packets" that are emit- 
ted stochastically both by the source and by the surrounding 
envelope. The optical depth determines the location at which 
the packets interact while their albedo defines the probability of 
either scattering or absorption. In the original scheme (scheme 
1) the source and the envelope photon packets are emitted sep- 
arately. At first the grains re-emit according to the absorbed 
source radiation. Then dust reemission takes also into account 
the envelope emission radiation field. Reemission by the dust is 
repeated as long as the difference between the input and the out- 
put energy is larger than a chosen threshold. However, the dust 
reemission, i.e. the repetition of the Monte Carlo experiment, 
is time consuming. An alternative possibility (scheme 2) is to 
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Fig. 1. Optical data for spherical astronomical silicate grains 
having a radius of 0. 12 fim (Draine & Lee 1984). Note that scat- 
tering dominates between 0.2 and 1 fim for this type of grains. 
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Fig. 2. Density structure perpendicular to the disk midplane 
and centered on the star for the model with t v = 1. Values 
are normalized to the maximum density. The contours provide 
0.10, 0.19, 0.28, 0.38, 0.48% of the maximum. 



store all radiation exchanges within the envelope. In this case 
the Monte Carlo experiment can be carried out once for all 4 , but 
a large amount of computer memory is needed. A drawback of 



4 This is only valid for opacities explicitely independent on the tem- 
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Table 2. Main features of the codes. 



Feature 


MC3D 


MCTRANSF 


RADICAL 


RADMC 


STEINRAY 


3D 


+ 








+ 


anisotropic scattering 


+ 


+ 




+ 


+ 


arbitrary grid geometry 


+ 


+ 








grain size distribution 


+ 


+ 


+ 


+ 


+ 


multiple dust species 


+ 




+ 


+ 


+ 


images 


+ 


+ 


+ 


+ 


+ 


polarization maps 


+ 










global error control 






+ 




+ 


multiple/extended heating sources 


+ 




+ 




+ 


dust evaporation included 




+ 


+ 






acceleration for high r (ALI/Ng/other) 


+ 


+ 


+ 




+ 


parallel version 




+ 


+ 







these two schemes is that the input luminosity is not automati- 
cally conserved during the simulation. This becomes a serious 
problem for configurations with very high optical depths which 
therefore usually need a larger number of iterations. A solution 
has been found by Bjorkman & Wood (2001 ) in the so-called 
immediate reemission (scheme 3). In this case only source pho- 
ton packets are emitted and followed in their interaction loca- 
tions. When a packet is absorbed, its energy is added to the en- 
velope and a new packet is emitted immediately at a frequency 
which takes into account the modified envelope temperature. 
This method does not require any iteration and implicitly con- 
serves the total energy. Another improvement of the standard 
MC procedure has been proposed by Lucy ( 1999 1 to treat ex- 
tremely optically thin configurations. This approach considers 
the absorption not only at the end points of the photon path but 
also in between. 

Grid-based codes solve the RT equation on a discrete spa- 
tial grid. The grid can be either determined during the sim- 
ulation or generated before starting the computation. The 2 
RT grid-based codes we compare, namely RADICAL and 
STEINRAY, use the second approach (see Sect. 13.2.31 and 
I3.2.5> . Among the schemes applied to solve the RT problem 
two are the most used: the so-called "Lambda Iteration" (see 
e.g. Collison & Fix fTWTl Efstafhiou & Rowan-Robinson fTWTl 
and the "Variable Eddington Tensor" (Mihalas & Mihalas 
[19841 Malbet & Berto utfT99ll Stone et a l.fT992l Kikuchi et al. 
120021 Dullemond et al. l200l Dullemond l2003l . The "Lambda 
Iteration" mode suffers from the same convergence problems 
as the standard MC which is based on scheme 1. Thus, the 
improved "Accelerated Lambda Iteration" (e.g. Rybicki & 
Hummer 1991 1 method is more widely applied. The "Variable 
Eddington Tensor" mode is more robust than the "Lambda 
Iteration" and usually converges faster. Moreover, it has been 
proven that it works properly even at extremely high optical 
depths. Integration of the formal RT equation can be done in 
a straightforward way by applying the "Long Characteristics" 
algorithm. This method is accurate but turns out to be costly 
in CPU time. A more efficient way is based on the "Short 



Characteristics" algorithm (Mihalas et al. 1978, Kunasz & 
Auer fT988l . 

Each of the solution algorithms has its advantages and 
drawbacks. In MC methods, a photon is propagated through 
the calculation domain and its scattering, absorption, and re- 
emission are tracked in detail. This allows to treat very com- 
plicated spatial distributions, arbitrary scattering functions and 
polarization. Drawback is the presence of a random noise in the 
results. This noise can be reduced by increasing the number 
of used photon packages and by including deterministic ele- 
ments in the MC experiment (Niccolini et al. l2003> . Grid-based 
solvers are less flexible than MC codes but have the advantage 
not to involve random noise. 

3.2. Code description 

In the following sections we briefly describe the RT codes par- 
ticipating in the 2D benchmark. A summary of their main fea- 
tures is provided in Table |2 

3.2.1. MC3D 

MC3D is a 3D continuum RT code. It is based on the MC 
method and solves the RT problem self-consistently. MC3D 
is designed for the simulation of dust temperatures in ar- 
bitrary dust/electron configurations and the resulting observ- 
ables: spectral energy distributions, wavelength-dependent im- 
ages and polarization maps. 

For the estimation of temperatures either the standard 
scheme 1 (see Wolf et al. I1999bl Wolf & Henning 2000 1 
or the immediate reemission concept (scheme 3) can be ap- 
plied. For this benchmark project, the scheme 3 is used to treat 
properly the more optically thick models. Optically very thin 
configurations, such as the atmosphere/envelope described in 
Sect. 12. 21 are easily computed by the method proposed by Lucy 
(1999 1. Furthermore, the efficiency of MC3D is increased by 
(a) the fast photon transfer and (b) wavelength range selec- 
tion concept (see Wolf & Henning 2000 for details), and (c) 
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the enforced scattering mechanism as described by Cashwell 
& Everett (IT9591 

Previous applications of MC3D cover feasibility studies of 
extrasolar planet detections (Wolf et al. 2002), the RT in the 
clumpy circumstellar environment of YSOs (Wolf et al. 1998 1, 
polarization studies of T Tauri stars (Wolf et al. 200jk)> AGN 
polarization models (Wolf & Henning 1999a I, a solution for 
the multiple scattering of polarized radiation by non-spherical 
grains (Wolf et al. I2002L and the inverse RT based on the 
MC method (Wolf l200Tb '). Executables of MC3D (V2) can be 
downloaded for several model geometries and platforms from: 
http://www.mpia-hd.mpg.de/FRINGE/SOFTWARE/mc3d/ 



(current US mirror page: http://spider.ipac.caltech.edU/staff/swolf/mc3d/i. 



3.2.2. MCTRANSF 

MCTRANSF solves multi-dimensional continuum RT prob- 
lems in dusty media by means of a MC method. It has been 
originally developed by Lopez et al. ( 1995 ). So far, the code has 
been used for the empirical modelling of several circumstellar 
envelopes of post AGB-stars of different types (e.g. Lopez & 
Perrin l2000t . including multi-scattering effects. 

Currently, only spherical symmetric (ID) and axisymmet- 
ric (2D) problems can be considered, but an extension to 3D 
is possible and straightforward. Several improvements of the 
standard MC procedure have been recently included (Niccolini 
et al. 12003b in order to avoid the usual increase of the noise 
level which typically occurs in extremely optically thin or opti- 
cally thick situations. The concept suggested by Lucy ( 1999 1 is 
implemented to treat very optically thin cases. Optically thick 
configurations are tackled by the inclusion of several determin- 
istic elements for the treatment of the absorption during the 
photon propagation phase, forcing the absorption to take place 
all along the rays. The temperature structure of the medium in 
radiative equilibrium is found by applying scheme 2. The con- 
vergence is found to be rapid, even in optically thick situations, 
but needs a large amount of computer memory, because the pri- 
mary MC information must be stored source-dependently. 

As a result of the combination of all these measures, 
MCTRANSF is capable to simultaneously model optically thin 
and optically thick parts of the model volume with about the 
same accuracy. Thus, both clumpy media and discontinuous 
opacity structures can be handled. MCTRANSF is able to ar- 
rive at numerical solutions for RT problems even in case of very 
large optical depths (e.g. for disk configurations). Parallelised 
versions of the code have been developed for a Cray T3E 1200 
and for systems supporting the OpenMP application program 
interface. All these versions use shared memory systems. 

3.2.3. RADICAL 

The core of the code RADICAL is a lambda operator subrou- 
tine based on the method of "Short Characteristics", imple- 
mented on a polar grid by Dullemond & Turolla (2000). Using 
this subroutine as the main driver, RADICAL offers two modes 
of operation: a simple Lambda Iteration mode and a Variable 
Eddington Tensor mode. In this paper we use the latter because 
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Fig. 3. Grid adopted by MCTRANSF to store the temperature 
resulting from the RT simulations. Similar spherical grids are 
also used by the other codes. 



of its faster convergence and capability of treating high optical 
depths. The Variable Eddington Tensor method is implemented 
in RADICAL as follows (Dullemond 2003 ). First, the primary 
stellar radiation field is propagated from the star outwards into 
the disk. Dust scattering is included in a MC fashion. The en- 
ergy absorbed by the disk in each grid-cell is then re-emitted 
as infrared (IR) radiation, which is treated as a separate radi- 
ation field. The 2-D transfer solution for this secondary radia- 
tion field is found by solving the frequency-integrated moment 
equations. The closure for these equations is based on the vari- 
able Eddington tensors and mean opacities computed with the 
Short Characteristics method of Dullemond & Turolla. At the 
end of the calculation a global check on flux conservation is 
made. For all the models discussed in this paper the error re- 
mained within 2%. 

3.2.4. RADMC 

RADMC is an MC code based on scheme 3. However, the 
original method of Bjorkman & Wood produces very noisy 
temperature profiles in regions of low optical depth, and 
requires a large number of photons (N ~ 10 7 ) for a smooth 
SED. These disadvantages have been solved in RADMC by 
treating absorption partly as a continuous process (Lucy 1999 
), and using the resulting smooth temperature profiles with a 
ray-tracing code to produce images and SEDs. These images 
and SEDs have a low noise level even for relatively few photon 
packages (N ~ 10 5 ). In addition, the frequency grid used for 
RADMC is not bound by the constraints set in the original 
method. This improved Bjorkman & Wood method works well 
at all optical depths, but may become slow in cases where the 
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Table 3. Resolution and number of photons for the different test cases. 



Code 


#r 


Ar 


#9 


A£> 


#Phot 


Test case 


Name 




[AU] 




[°] 


[xlO 6 ] 


T 550nra 


MC3D 


55 


0.03-141 


121 


1.5 


244 


0.1,1,10 


MC3D 


10 3 


0.07-4.1 


121 


1.5 


244 


100 


MCTRANSF 


48 


0.17-125 


40 


4.5 


1000 


0.1 


MCTRANSF 


48 


0.17-125 


40 


4.5 


800 


1 


MCTRANSF 


46 


0.18-130 


46 


2.8-5.3 


1000 


10 


MCTRANSF 


46 


0.18-130 


46 


2.8-5.3 


500 


100 


RADICAL 


60 


0.03-116 


62 


1.6-8.3 




0.1-100 


RADMC 


60 


0.03-116 


62 


1.6-8.3 


10 


0.1-100 


STEINRAY 


61 


0.12-109 


61 


1.3 




0.1-100 



optical depth is very large (t v about 1000). For the test cases in 
this paper the optical depths are low enough that this problem 
does not play a role. For more information on the code, see 
http : //www . mpa-garching . mpg . de/PUBLICATIONS/DATA/ 
radtrans/radmc/. 

3.2.5. STEINRAY 

STEINRAY is a grid-based code which solves the full 3D con- 
tinuum RT problem. A combination of ray-tracing and finite 
differencing of 2nd order on adaptive multi-frequency pho- 
ton transport grids is applied. Steinacker et al. (2002ai have 
shown that the use of 3rd order finite differencing is too time- 
consuming for 3D RT, while 1st order schemes introduce an 
unacceptable degree of numerical diffusion to the solution. 

The spatial grids are generated using an algorithm de- 
scribed in Steinacker et al. ( 2002c ). They are adaptive and opti- 
mized to minimize the 1 st order discretization error hence guar- 
anteeing global error control for solutions of radiative transfer 
problems on the grid. Since the use of one single grid for all fre- 
quencies leads to large discretization errors, STEINRAY calcu- 
lates individual grids for each frequency to use the global error 
control of the grid generation method. Minimization of the grid 
point number is possible in regions where the optical depth be- 
comes large allowing for treatment of applications with opti- 
cal depth of any value. Contrary to former treatments, the full 
frequency-dependent problem is solved without any flux ap- 
proximation and for arbitrary scattering properties of the dust. 
For the direction discretization, equally spaced nodes on the 
unit sphere are used along with corresponding weights for the 
integration derived by evaluating special Gegenbauer polyno- 
mials in Steinacker et al. ( 1199 6). The temperature distribution 
is calculated by an Accelerated Lambda Iteration between the 
radiative transfer equation and the local balance equation. The 
program is designed to provide spatially resolved images and 
spectra of complex 3D dust distributions and allows for multi- 
ple internal and external sources (Steinacker et al. l2003> . 

Recently, a 2D version of the program has been developed 
and was used for this benchmark. The grids are similar to the 
spherical grids shown in Fig. [5] The 2D version uses ray-tracing 
to solve for the intensity in all directions. 



3.3. Reliability of the codes in 1D 

All the RT codes participating in the benchmark have been 
already tested in ID spherically symmetric configurations. 
Results from MC3D have been compared with those calculated 
with the RT code of Chini et al. ( 1986 1 and with the code of 
Menshchikov & Henning ( 1997 1. In both cases differences be- 
low 1 % even in the case of high optical depth have been found 
(Wol f et al . I1999bi . MCTRANSF has been tested (Niccolini 
et al. 2003 1 against the ID code written by Thibaut le Bertre 
and based on the work of Leung ( 1976). For the case of optical 
depth equal to 10 at 1 fim (with k oc i, and isotropic scat- 
tering) the agreement between MCTRANSF and le Bertre's 
code is better than 1%. In the case of RADICAL, tests have 
been performed by comparing the results with those from 
TRANSPHERE 5 , a 1-D variable eddington factor code tested 
against the similar "code 1" of Ivezic et al. J1997i . Steinacker 
et al. (2003) used the ID benchmark provided by Ivezic et al. 
( 119971 ) to test the 3D version of STEINRAY. Agreement be- 
low 1 % has been found both for the emerging temperature and 
SEDs. 

3.4. Details on the 2D RT computation 

The three MC codes involved in the 2D benchmark, namely 
MC3D (Sect. l3~2~Tt and MCTRANSF (Sect. l3~Z2t and 
RADMC (Sect. l3~2~4l ). choose a spherical grid to store the tem- 
perature resulting from the RT simulations. Radially the steps 
are logarithmic in order to properly resolve the innermost dense 
region of the disk. The number of radial points is kept be- 
low 100 for all the codes but MC3D, for which we tried a 
two weeks long computation for the most optically thick case 
(see Table [3}. Differences in the SED between this computa- 
tion and the one with 55 radial points and lower number of 
photons are discussed in Sect. 14.41 MC3D adopts an equally 
spaced grid in vertical direction (1.5° resolution, see Table|3j, 
while MCTRANSF and RADMC choose a resolution decreas- 
ing with the distance from the disk midplane (see Fig. [3j. 
Similar grid geometries are also used by the two grid-based 
codes RADICAL and STEINRAY. RADICAL (Sect. l3~2~5l 

5 http://www.mpa-garching.mpg.de/PUBLICATIONS/DATA/ 
radtrans/ 
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makes use of the same grid as RADMC while STEINRAY 
(Sect. 13.231 has comparable resolution in vertical direction but 
larger cells in the inner disk (0.12 AU, see Table[3}. 

We note that the resolution given in Table [3J together with 
the number of photons for the MC codes, are at the limit of 
the computing capabilities for most of the codes. MC3D needs 
about 1 Gby memory. The temperature resulting from these 
test cases is obtained in 1-2 days. Computing the SED requires 
about a week for all the models, but for the most optically thick 
one for which we try a longer computation. MCTRANSF needs 
a large amount of computer memory to store all radiation ex- 
changes. The necessary memory goes as the square of the num- 
ber of cells: 46x46 cells is actually the technical limit (~ 4Gby) 
on our 4-processors (ev67 at 1 GHz) HP-compaq ES45 server. 
The runtime for the most optically thick case is about 2 weeks. 
Results from RADICAL and RADMC are obtained in less than 
a day for all the test cases. However, the actual spatial resolu- 
tion of RADICAL cannot be doubled due to not sufficient com- 
puter memory. To produce the final spectrum RADMC uses the 
ray-tracing module from RADICAL, and for that reason is also 
limited to the maximum spatial resolution that can be achieved 
by RADICAL. This is a mere technical problem. In the 2D 
version of STEINRAY, the code uses about 1 Gby memory. To 
obtain convergence in the temperature iteration of 1 % for the 
case t v = 100, the code runtime is about a week. 

4. Results 

4.1. Approximate solution for optically thin 
configurations 

In the case of configurations being optically thin for all rel- 
evant wavelengths, heating of the dust particles is dominated 
by the stellar radiation. When re-emission of the dust particles 
can be neglected and scattering is only included as extinction 
term, the dust temperature can be easily determined from eq. 
without any coupling to eq. Q. We use this approximate 
semi-analytical solution to test independently our RT compu- 
tations for the most optically thin case and to check the cor- 
rectness of the density setup in the other more optically thick 
models (see Sect.lOandPDl. 

The assumptions we discussed above simplify eq. in the 
following way 
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with T&(R, 9) being the disk temperature at the location (R, 9) 6 
and Q (A) the absorption efficiency factor. We perform the 
integration at the same wavelengths as those adopted by the RT 
codes (see Sect. 12. 21 . The term B A (T t ,R, 9) represents the black 
body emission from the star corrected for the extinction 



B e A (T„,R,9) = B A (T*)e 



-Tra^Gf +G. s r) J? dR'p(R',6) 



(6) 



where a is the dust radius, <2 sca (/l) is the scattering efficiency 
factor and p(R', 9) is the density distribution given in eq. @ but 



6 R is the distance from the central star in spherical coordinates and 
6 is the polar angle as measured from the disk midplane 



Diff. [%] 
5 








i=12.5° " 







VK 






-5 











0.1 



1.0 



10.0 



100.0 



1000.0 



10000.0 



X [um] 



Fig. 4. Differences of the most optically thin model from the 
semi-analytical solution (see Sect. 14. II . Upper panel: differ- 
ences in radial temperature for an angle 9 near to the disk 
midplane. Lower panel: differences in the SED for an almost 
face-on disk (disk inclination equal to 12.5°). For both panels, 
solid lines give the difference of MC3D, dot-dashed lines of 
MCTRANSF, dashed-dot-dot-dot of RADICAL, dotted lines 
of RADMC, and dashed lines of STEINRAY from the semi- 
analytical solution. 

here expressed in spherical coordinates. The argument of the 
exponent represents the optical depth t a (R, 9) at the distance R, 
9 from the central star. Once the optical depth is determined, 
the extincted black body emission can be substituted in eq.Q 
and the dust temperature can be easily computed. To obtain the 
flux density F A at a distance equal to the stellar radius we need 
to integrate the power emitted by each grain over the entire 
volume. If we express the power emitted by one grain as 



Pl(R,9) = 4na 2 Qf s B A [T d (R,9)] 



(7) 



the flux density can be obtained by solving the following inte- 
grals 



2tt 
4nR: 



/"•Rom f*n 

I I ' 



d9'dR'P s A (R',9')p(R',9')R' 2 cos(9') (8) 



here the factor 27r comes from the integration in <p which is the 
azimuthal angle in the x — y plane. 

A first rough estimate of the correctness of this approach 
can be done by evaluating how much the disk temperature T& 
would increase because of secondary emission re-absorption 
events. The grain emissivity is proportional to T% and, in case 
of optically thin emission, to the optical depth. Our most op- 
tically thin model has a maximum IR optical depth of 0.01 at 
10 fim in the disk midplane, where most of the dust is confined. 
If an emitted IR photon were re-absorbed by the disk, its tem- 
perature would increase by the quantity (1 + Tir) 25 . Thus, the 
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Fig. 5. Radial temperature (upper panel) and percentage of dif- 
ference (lower panel) among the codes for the most optically 
thick model. RADICAL is taken as reference code. The radial 
cut is made for an angle 6 near to the disk midplane. Diamonds 
give the radial dependence in case of long wavelengths and op- 
tically thin emission. In the upper panel, solid lines are the re- 
sults fromMC3D, dot-dashed lines from MCTRANSF, dashed- 
dot-dot-dot from RADICAL, dotted lines from RADMC and 
dashed lines from STEINRAY. In the lower panel, solid lines 
give the difference of MC3D, dot-dashed lines of MCTRANSF, 
dotted lines of RADMC and dashed lines of STEINRAY from 
RADICAL. 



temperature obtained by the semi-analytical approach can be 
considered correct within 0.26% in a first approximation. The 
corresponding emergent flux has an uncertainty which is about 
four times larger (about 1%). However looking at the SED dif- 
ferences (see Fig.@}, it is clear that deviations due to scattering 
(treated correctly in the numerical codes) are important too. 
A more realistic estimate of the error on the semi-analytical 
approach should indeed consider the effect of scattering apart 
from damping the stellar flux. 

4.2. Temperature 

All the codes correctly reproduce the shape of the temperature 
distribution, both its radial and vertical dependence. 

In oder to test independently the most optically thin case 
(highest optical depth in the disk midplane of t v = 0.1), we 
use the semi-analytical solution derived in Sect. 14. II The upper 
panel of Fig. @] shows the percentage of difference 7 between 
the semi-analytical solution and any other code. Radial cuts 
are plotted for an angle 6 near to the disk midplane where the 



with difference we mean: 



(g-b) 
b ' 



Here b stands for the reference 



solution/code while a for any other code. 



Fig. 6. Vertical temperature and percentage of difference 
among the codes for the most optically thick model and for 
a distance r in the midplane equal to 2 AU from the central 
star. RADICAL is taken as reference code. The nomenclature 
is the same as for Fig. |5] except for the diamonds which give 
the temperature behaviour following the semi-analytical ap- 
proach described in Sect. 14.11 Note that all the RT codes have 
the same turnover in the temperature distribution at the location 
predicted by the semi-analytical solution. 



influence of scattering is the largest. The differences between 
the semi-analytical solution and the RT codes is always smaller 
than 1%. 

Temperature distributions for the models with r v =1 and 
10 agree better than 10% for all the cases we examine. The 
disk model with the highest optical depth is the most difficult 
to treat for the RT codes. In Fig. [5] and [6] we show radial and 
vertical cuts at the disk locations where deviations from the 
codes are expected to be higher, i.e. near to the disk midplane 
and close to the central star. The radial temperature is plotted 
for an angle 9 equal to 2.5° from the disk midplane while the 
vertical temperature is given for a distance equal to 2 AU from 
the central star. In the upper panel of Fig. [5] we also superim- 
pose the temperature dependence for the optically thin regime 
at long wavelengths. In this regime the temperature distribu- 
tion depends only on the dust properties and can be approxi- 
mated by T(r) oc r' 2 ^ + ® (Evans 19941. Here /3 corresponds to 
the index of the dust absorption coefficient at low frequencies 
(a:^ 1 oc v 13 ). For Draine & Lee silicates /3 is equal to 2, leading 
to an exponent of -0.33 in the temperature relation. 

The upper panel of Fig. [6] provides (in diamonds) the ver- 
tical temperature profile from the semi-analytical approach 
described in Sect. 14.11 The semi-analytical solution has the 
turnover point from optically thick to optically thin (the place 
where the temperature suddenly starts to drop) around 19° 
from the midplane. Since the solution provided by the RT 
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Fig. 7. SED for two disk inclinations i as given on top of each panel. Each curve provides the mean value of the five RT simulations 
for the four computed models with different optical depth. The midplane optical depth is given in parenthesis labeling each curve. 
In both panels solid lines show results for the most optically thin disk, dotted lines for a disk having r v = 1, dot-dashed lines 
for a disk with r v = 10 and dashed lines for the most optically thick model. Diamonds provide the black-body emission from 
the naked star. The slope of the SED at long wavelengths depends only on the dust properties and is plotted in each panel with a 
solid line. 



codes should have the same behaviour, we used the semi- 
analytical approximation to check the correctness of the density 
setup. The disk midplane calculated with the semi-analytical 
approach is naturally cooler than the real disk because the ap- 
proximation neglects heating of the disk from dust re-emission. 
The outer regions of the real disks are also warmer because of 
those photons scattered far from the midplane. At a distance r 
of 2 AU from the star and exactly in the midplane the real tem- 
perature is a factor of about 1.3 higher than that given by the 
semi-analytical solution. 

The lower panels of Fig. [5] and |6] provide the percentage 
of difference among the codes taking RADICAL as reference 
code. The radial temperatures agree better than 5% in most 
of the disk, going from 1.2 AU to 200 AU. Around 1.1 AU 
MCTRANSF deviates slightly more than 10%. STEINRAY 
shows 10% deviations at the inner boundary and slightly higher 
deviations (but always less than 15%) far from the star, at about 
900 AU. The vertical cut at 2 AU shows an agreement better 
than 2.5% till 10° from the disk midplane. Closer to the disk 
midplane, deviations are larger for MC3D and STEINRAY but 
always smaller than 4%. 

4.3. Spectral Energy Distribution 

The emerging spectra for the four models having different opti- 
cal depths are shown in Fig.[7]at two disk inclinations. The left 
panel provides the results for an almost face-on disk (disk incli- 
nation i equal to 12.5°) while the right panel gives the result for 



an almost edge-on disk (i = 77.5°). Each curve represents the 
mean value of the five RT simulations for the specific model, 
whose optical depth is given in parenthesis above the curve. On 
the y axis, we plot A in [W m~ 2 ] where F,\ is the flux den- 
sity at a distance equal to the star radius. We also superimpose 
in diamonds the black body radiation arising from the star in 
order to visualize how efficiently the circumstellar disk repro- 
cesses the stellar energy. We note that all the codes have the 
correct slope at long wavelengths. This slope depends only on 
the dust properties and is plotted as solid line in both panels 
(AFj oc A~ 5 ). At 0.55 fim the drop in luminosity amounts to 
about a factor of 20 going from the most optically thin to the 
most optically thick model and for a disk inclination of 77.5°. 

Since the differences among the codes are too small to be 
visible in a logarithmic plot, we provide separately the percent- 
age of difference for the four models and for three disk incli- 
nations (see Fig. [8). RADICAL has been chosen as reference 
code. For the most optically thin case, we also compare our 
results with the semi-analytical approach (see Fig.0}. We find 
that the agreement of the codes with the semi-analytical solu- 
tion is always better than 8%, with the largest deviations around 
0.3 and 40 fim. In the range 0.2-0.7 jum all the codes predict 
higher flux in comparison to the semi-analytical solution, while 
between 10-200 fim a lower flux is obtained. These deviations 
arise because the semi-analytical approach includes scattering 
only as an extinction term. From the numerical RT calculations 
it is clear that some photons are scattered thus enhancing the 
flux between 0.2 and 0.7 /mi. We note that this wavelength 
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Fig. 8. Percentage of difference in the SED between the codes. RADICAL is taken as reference code. Solid lines give the dif- 
ference between MC3D and RADICAL, dot-dashed lines between MCTRANSF and RADICAL, dotted lines between RADMC 
and RADICAL and dashed lines between STEINRAY and RADICAL. 



range is exactly where small astronomical silicate grains have 
the largest scattering efficiency (see Fig.0. Therefore, devia- 
tions peaking at 0.3 yt/m are simply explained by the particular 
optical data chosen for this benchmark. Those photons which 
are scattered cannot contribute to heat the disk. This explains 
why RT codes predict a lack of emission at longer wavelengths. 
To understand why the largest deficit of photons is around 40- 
50 jt/m, we first compute the temperature at which most of 
the disk mass emits (mass average temperature) and then the 
corresponding wavelength. For the wavelength calculation we 



need to take into account the grain emissivity (Evans 1994 1. 
We find a mass average temperature of 40 K which translates 
into a wavelength of 50 fim at the maximum emission. This 
wavelength is well in agreement with the deviations shown in 
Fig-El For comparison, the RT codes agree better than 1.5% 
at wavelengths shorter than 10 //m for this particular test case 
(6 = 12.5°). At longer wavelengths the results show a bit more 
scatter but the agreement is always better than 3%. In Fig. [H] 
first panel, we also show the percentage of difference for two 
other disk inclinations, namely 42.5° and 77.5°. In both cases 
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the agreement is better than 2% at all wavelengths for all the 
codes but MCTRANSF, for which slightly higher deviations 
(about 2.5%) are present at longer wavelengths. 

As the optical depth in the midplane increases, the RT prob- 
lem becomes more difficult to solve. Because of the chosen disk 
geometry, most of the disk mass is located near to the mid- 
plane and close to the disk inner boundary. Our comparison 
shows that agreement among the codes is always better for an 
almost face-on case and 42.5° disk inclination (first two pan- 
els of Fig.[S}, than for an almost edge-on disk (lower panels of 
Fig-EJ- For the models with r v =1 deviations among the codes 
are smaller than 9%. For the disk with midplane optical depth 
of 10 and 100 and inclinations of 12.5 and 42.5°, differences 
do not exceed 10%. For the almost edge-on configurations the 
most difficult regions to treat are those where scattering dom- 
inates and at wavelengths around 10 fim. In the IR, opacities 
change strongly and the modified Planck emission peaks in the 
inner disk regions (between 1 and 2 AU). Therefore, the nu- 
merical simulations are particularly sensible to the resolution 
of the inner parts. Deviations in the IR are partly due to the 
different resolution adopted by the codes (see also Sect. I4.4i . 
Scatter at visible and near-infrared wavelengths for MC3D, 
MCTRANSF, and RADMC is simply statistical noise, typical 
of MC simulations. This scatter becomes more prominent at 
high optical depths. We also note that MC3D and RADMC 
have the same trend for wavelengths larger than 10 //m and 
the model with r v =100: they both estimate a larger IR emis- 
sion than RADICAL with peaks at -70 /urn for MC3D and at 
~ 100 /im for RADMC. A strong deviation from the other codes 
is shown by STEINRAY just after the 9.8 /mi silicate feature. 
However, one should note that apart from the discussed features 
the overall agreement of the SEDs is better than 10% for all the 
codes even for the almost edge-on disk and the most optically 
thick test case. 



4.4. Tests for various spatial and frequency resolutions 

We used the MC code MC3D to test the dependence of our 
results on the grid adopted to store the emerging temperature. 
Since deviations due to different temperature sampling are ex- 
pected to be larger for more optically thick configurations, we 
investigate our most optically thick test case t v = 100. Different 
grids, covering radial resolutions from 2.7 AU up to 0.03 AU in 
the inner disk and vertical resolutions from 1.5° to 5.8°, have 
been inspected. In Fig. [9] we report results for five relevant 
cases. The number of radial and vertical subdivisions (#r and 
#8), as well as the resolution (Ar and A6) and total number of 
photons (#Phot) for these cases are provided in Table 0] The 
RT equation is solved for 61 wavelengths, the same assumed 
in all the previous simulations. The number of photons is set to 
4xl0 5 per wavelengths to limit the runtime to 2 days on a PC 
with 4 Gby memory, 2.4 GHz clock. Only in one case (mod4) 
we let the computation run for about two weeks in order to re- 
duce as much as possible the photon noise at short wavelengths. 
The comparison shows that the vertical spacing does not influ- 
ence too much the results: we report deviations smaller than 
1.5% for the three disk inclinations between the models with 



Table 4. Relevant models for the spatial resolution tests 
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Fig. 9. Spatial resolution tests using the MC code MC3D. On 
the y-axis we plot the percentage of difference between the 
emerging SED of modO and any other model in Table |4] Solid 
line: difference between modO and modi. Dot-dot dahsed line: 
difference between modO and mod2. Dotted line: difference be- 
tween modO and mod3. Dashed line: difference between modO 
and mod4. 

vertical resolution 1.5° and 5.8° (solid line in Fig.|9j- On the 
other side, changes in the radial grid strongly affects the IR 
emission: when going from modO to mod2 results still do not 
differ more than 5% but the coarse inner grid of mod3 causes 
deviations larger than 20%. Model mod4 is the state-of-the-art 
for our computer capabilities. The grid has a good resolution 
also in the outer region of the disk and we use 10 times more 
photons per wavelengths. The percentage of difference between 
this model and modO in the IR regime amounts to less than 6%. 
Deviations of about 10% in the optical range are due to differ- 
ences in the number of photons. 

To test the influence of the frequency resolution on our re- 
sults we run the 2 MC codes MC3D and RADMC doubling 
the number of walenghts where to solve the RT equation. The 
other three codes could not take part to the comparison because 
of not enough computer memory. To have reasonable runtime 
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Fig. 10. Frequency resolution test using the MC codes MC3D 
and RADMC. On the y-axis we plot the percentage of differ- 
ence between the emerging SED of the codes for two models, 
one sampling the frequency space with 61 (dotted line) and the 
other with 122 (dashed line) logarithmically distributed points. 

for MC3D (a couple of days), we restrict ourselves to the case 
r v = 10 and we lower the number of photons in comparison 
to Table (half in the case of RADMC and 4 times lower in 
the case of MC3D). The lower number of photons introduces 
larger scatter at short wavelengths. Fig.llOlshows the percent- 
age of difference between MC3D and RADMC for the model 
with 61 wavelenghts (dotted line) and the new model with 122 
wavelenghts (dashed line). Apart from the expected larger de- 
viations at short wavelengths, the agreement between the two 
codes improves not more than 2% around 10//m. Thus, we con- 
clude that the nature of the IR deviations plotted in Fig. [8] is not 
due to coarse frequency sampling but to the different grid res- 
olutions adopted by the codes (especially in radial direction) 
together with cumulative numerical errors. 

5. Discussion and Conclusions 

Before presenting our findings, we briefly discuss the general 
features of the computed SEDs for different viewing angles (i) 
and optical depths (r). As mentioned in Sect. 12.21 optical depths 
are measured in the disk midplane from the inner to the outer 
boundary. Thus, the optical depths we refer to are the highest 
in the disk. 

Both panels of Fig.Qshow clearly that the far-infrared re- 
gion (long ward 100 /vm) remains unaffected when varying the 
viewing angle. On the other hand, the short-wavelength part of 
the spectrum is strongly modified. When the disk is seen face- 
on the spectrum is dominated by unattenuated stellar radiation. 
As the disk inclination increases, more and more of the stellar 
flux is extincted by the dust in the disk. For the t v = 100 case 



at i = 77.5° this reduction amounts to a factor of e in the 
visual. However, due to the high albedo of the dust grains, a 
large fraction of the stellar radiation is scattered above the disk 
into the line of sight. Dust scattering is also responsible for the 
excess of emission at stellar wavelengths seen at small disk in- 
clinations (see left panel of Fig. especially for high optical 
depths). The optical depth also affects the strong 10 //m feature 
produced by the SI - O stretching. While in most cases the fea- 
ture appears strongly in emission, for the t v = 100 test case 
at i = 77.5° the feature appears in absorption (see Fig. fright 
panel). The 20 fim feature is much weaker than the 10 /im, but 
it is already visible for the model with t v = 1 . All these features 
are in agreement with earlier radiative transfer computations of 
disks (e.g. Efstathiou & Rowan-Robinson 1990 Menshchikov 
& Henning H997l . 

Our aim in this paper is to provide benchmark solutions for 
the 2-D continuum radiative transfer problem in circumstellar 
disks. The problems we present have optical depths up to 100, 
which is actually the limit of current computational capabili- 
ties for most of the codes. The corresponding total dust mass 
in the disk of about 0.01 solar masses covers most of the ob- 
served disks around low mass stars. For more massive disks 
around intermediate and high mass stars as well as tori obscur- 
ing active galactic nuclei, the numerical strategies have to be 
modified, using e.g. the diffusion approximation for high opti- 
cal depths. We used five independent radiative transfer codes 
that implement different numerical schemes. We compare both 
the resulting temperature structure and the emerging SEDs. For 
the lowest optical depth case (t v = 0.1) we also compared the 
results against a semi-analytic solution which treat scattering 
only as extinction term. The other three cases (r v = 1,10, 100) 
cannot be solved in a semi-analytic way, since multiple scatter- 
ing and absorption-reemission events strongly affect the solu- 
tion. 

We find that the overall shape of the temperature distribu- 
tion and of the emerging SEDs is well reproduced by all the 
codes. Differences in the temperature are smaller than 1 % for 
all the codes in the most optically thin case. Even for the most 
optically thick model, differences in the temperature remain 
below 15%. As for the SEDs, deviations among the codes are 
smaller than 3% at all wavelengths and disk inclinations for the 
most optically thin model. For the models with t v = 1 and 10 
at all disk inclinations and for the most optically thick case for 
disk inclinations of 12.5 and 42.5°, differences do not exceed 
10%. Only for the most optically thick case and an almost edge- 
on disk differences around 10 fim exceed 20% in the case of 
STEINRAY. We stress that this is the case for which the numer- 
ics is the most difficult: the codes have to treat both a very op- 
tically thin atmosphere and a thick disk midplane. Independent 
tests using two of the MC codes show that the frequency res- 
olution cannot account for the infrared deviations among the 
codes in the almost edge-on disk and the most optically thick 
model. Grid resolution especially in radial direction together 
with cumulative numerical errors play a major rule. The pre- 
sented results provide a robust way to test other continuum RT 
codes and demonstrate the possibilities of the current computa- 
tional capabilities. Temperature distributions and SEDs for all 
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the test cases are available at the web site: 
http://www.mpia.de/PSF/PSFpages/RT/benchmark.html 
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